r*D-R137  628 

|  UNCLASSIFIED 


FINITE  ELEMENT  ANAL VS IS  FOR  COHESIVE  SOIL  STRESS  AND  1/1 

CONSOLIDATION  PROBLE .  .  OJ)  CALIFORNIA  UNIV  DAVIS 
L  R  HERRMANN  ET  AL.  DEC  83  NCEL-CR-84.  006 
N62582-83-M-T062 


F/G  8/13 


NL 


J*Vv*5‘  * 

V- 


cr  84.00^;?  wm 

NAVAL  CIVIL  ENGINEERING  LABORAtc||^ 
Port  Hueneme,  California  > 


Sponsored  by  ’  v } 

NAVAL  FACILITIES  ENGINEERING  COMMAND 


FINITE  ELEMENT  ANALYSIS  FOR  COHESIVE  SOIL,  STRESS 
AND  CONSOLIDATION  PROBLEMS  USING  BOUNDING 
SURFACE  PLASTICITY  THEORY 


December  1983 


An  Investigation  Conducted  by 
UNIVERSITY  OF  CALIFORNIA,  DAVIS 


r  \ 

j  \ 

.  *rr.  r  o 


V* 


<&$  .  ’ 


N625SH1J-M-T062 


Approved  for  public  release;  distribution  unlimited. 


DTIC  FILE  COPY 


84  0 


047 


$ 


<a«U  TO  i«  lum  * 


111 

Hill  lift 


St!  hid!  it 
] 

l,r  ll®r^5 

I  \  $$ 


i  i  ii , 
II  iiiil!  I 


i U 


If  •!  'hill  .1.  ill _ V« 


||I!ll  Still  til  illlill!! j|l 

jj  §?«*!  UsSsas  |  ,3:  |.„3i!slS  1 1| 


[jlii.ssIsSsssl.Sa 
i  till  „j| 

i  iiii  msi  sill 


K  i«u  v«Ti  s«  x- 


%*  Tljjl  MM  mtrnrn  0mm  M— e« 


|  UWT  DOCUMENTATION  PAGE 

|  CR  84.006 

mm 

4.  nTLC  8M04NKJ 

Finite  Element  Analysis  for  Cohesive 
Soil.  Stress  and  Consolidation  Problem! 
Using  Bounding  Surface  Plasticity 
Theory 

1.  Tv.,  0.  »«»0«T  ,  «|*M  COVI.KO 

Flpa). 

Jan  1983  -  Oct  1983 

T.  auTi4w« 

Leonard  R.  Herrmann 

Kyran  D.  Mlsh 

•.  GOKT..CT  OK  WUT  MUMOtHTO 

N62583-83-M-T062 

1.  UWt—m  o— —nation  «•■  uto  tMMiu 

University  of  California,  Davis 

-  SKraKiiJSSRP® 

Tr023.03m0l.002 

it.  CQMT<iOi.kiM«  orrtcc  mmm  a«o  ammii 

Naval  Civil  Bnglnaarlng  Laboratory 

»t.  *«»0*T  DAT* 

December  1983 

|  Port  Hu enema,  CA  93043 

VI  AHMCV  MAMK  •  AOOHtlrO  hmm  C**~M*4  0«(re) 

Ravel  Facilities  Engineering  Command 

200  Stovall  Street 

l»-  IICUKTV  ClUI.  (M  l*M| 

Unclassified 

Alexandria,  VA  22332 

••  OUTMiaUTIOM  IT.TCMNf  t*4  Mm 

Approved  for  public  release;  distribution  unlimited. 

IT  OlTTKlOuTlOK  TT.Tfut.T  (ml  mm  mrnmlwmtt  mmimtmrn  M  Wlmmk  H.  ,1  ttUmmwmt  M.  IM 

i*  luwiMa'Mi  MTU 

ff  a(v  SO*OI  fCwMSsl  am  re  verse  if  e  eras  eery  m4  »4en#ifr  be  bier*  mm*~) 

Finite  element,  computer  program,  geotechnical  engineering, 
soil  constitutive  lav 

M  IMTUCT  (dmmtmmm  mm  rrnmmmmm  mum  M  mmtmmmmr  mm*  i ow.fr  »r  mil  «M«) 

vThe  equations  governing  the  consolidation,  and  the  stress  and 
strains  states  for  soil  structures  are  reviewed  and  their 
historical  development  is  discussed.  Numerical  analysis  con¬ 
cepts  are  used  to  express  these  equations  in  incremental  fora. 

A  variational  statement  of  these  incremental  equations  is 
formulated  and  used  in  the  development  of  a  comprehensive 
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-^finite  element  analysis.  The  concepts  used  in  developing  the 
variational  statement  are  somewhat  different  from  those  used  by 
most  other  investigators  and  appear  to  offer  certain  advantages 
for  inelastic  formulations.  Finally  results  obtained  with  the 
finite  element  analysis  are  compared  to  known  solutions  with 
good  results. 

For  the  convenience  of  the  reader;,  the  total  report  on  the 
project  is  presented  in  four  parts.  As  noted  above  a  description 
of  the  consolidation  theory  and  certain  theoretical  features  of 
the  finite  element  analysis  are  described  in  the  body  of  the 
main  report  (CR  84.006).  The  second  part  (CR  84.007)  describes  j 
the  numerical  evaluation  o/  the  incremental  form  of  the  bounding  : 
surface  model.  Finally^ "user's  manuals"  for  the  2-D  and  3-D  j 
finite  element  programs  are  given  in  two  additional  reports 
(CR  84.008  and  CR  84.009).  I 
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1.  SCOPE  OF  PROJECT 

The  goal  of  the  project  is  to  develop  two  special  purpose  finite  element 
codes  for  the  analysis  of  cohesive-soil,  stress  and  deformation  problems  including 
consolidation  effects.  Specifically  the  codes  are  to  make  use  of  the  new 
comprehensive,  bounding  surface  plasticity  constitutive  model  for  cohesive  soils 
[1,2,3]. 

2.  INTRODUCTION 

The  analysis  is  limited  to  small  deformations  and  displacements,  and 
classical  consolidation  theory.  In  addition,  it  is  restricted  to  two  ideal  conditions 
of  saturation.  The  first  is  when  the  soil  is  completely  saturated  and  consideration 
is  given  to  the  development  and  dissipation  (due  to  water  flow)  of  excess  pore 
water  pressure;  included  in  this  case  are  ideal  undrained  conditions.  The  second 
case  assumes  that  the  degree  of  saturation  is  sufficiently  small  (or  ideal  drained 
conditions  exist)  that  the  pore  water  pressure  is  zero,  and  the  presence  of  water 
can  be  completely  accounted  for  by  the  increased  unit  weight  of  the  soil. 

Consolidation  theory  for  saturated  soils  is  well  established  and  can  be 
found  in  a  number  of  references  [e.g.,  4-8],  The  form  of  the  theory  used  in 
this  work  is  taken  from  [8]  with  only  slight  modification  and  some  changes  in 
notation;  the  theory  is  summarized  in  a  later  section. 

A  number  of  finite  element  analyses  have  been  developed  [6-12]  for  soil 
consolidation  problems;  most  have  been  limited  to  linear  elastic  material  behavior 
which  is  unrealistic  for  cohesive  soils  and  none  have  used  the  newly  developed 
bounding  surface  plasticity  theory.  Because  the  finite  element  concepts  employed 
in  the  programs  are  standard,  the  section  describing  the  analysis  will  be  brief 
in  nature. 
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3. 


A  BRIEF  LITERATURE  REVIEW  OF  THE  THEORY  AND  FINITE  ELEMENT 
APPLICATIONS  FOR  CONSOLIDATION 


This  section  reviews  the  theory  of  consolidation  and  examines  several 
representative  Finite  Element  analyses  for  its  evaluation.  It  is  not  an  exhaustive 
catalogue  of  research  in  consolidation  theory  and  analysis,  but  instead  attempts 
to  demonstrate  some  of  the  advantages  and  the  difficulties  ptesent  in  Finite 
Element  models. 

3.1  Consolidation  Theory 

The  birth  of  Soil  Mechanics  as  a  modern  engineering  discipline  occurred 
in  the  1920's,  when  the  Austrian  engineer  Karl  icrzaghi  proposed  his  theory  of 
the  consolidation  of  saturated  fine-grained  soils  under  applied  loads  [13]. 
Terzaghi  had  been  studying  the  phenomenon  of  the  reduction  in  void  space  of 
soils  underlying  foundations.  He  correctly  perceived  that  the  time-dependent 
settlement  from  consolidation  of  these  soils  was  due  to  the  flow  of  water  out 
of  the  soil  skeleton  as  the  voids  decreased  in  size.  The  permeability  of  the 
soil  dictates  the  rate  at  which  these  movements  take  place.  The  soil  skeleton 
thus  acts  like  a  large  sponge  in  response  to  an  applied  load. 

Many  of  the  most  important  features  of  com  fil  iation  can  be  motivated 
by  considering  a  greatly  simplified  mechanical  model,  the  'spring  analogy',  Fig.  1. 

Consider  a  cylinder  which  contains  a  piston,  valve,  and  elastic  spring, 
Fig.  la.  This  cylinder  is  filled  with  an  incompressible  fluid.  If  a  force  is 
applied  to  the  piston  with  the  valve  open,  the  force  is  initially  carried  by  the 
fluid,  Fig.  lb.  With  time,  however,  the  fin  d  drains  from  the  cylinder  under 
the  applied  force,  and  more  of  the  force  is  carried  by  the  spring,  Fig.  Ic. 
Finally,  a  new  equilibrium  position  is  reached,  Fig.  Id,  where  all  of  the  force 
is  carried  by  the  spring,  and  the  excess  fluid  pressuie  drops  to  zero. 
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Fifuri  1.  Sprint  A«loiy  For  Cousolidatia 


This  simple  analogy  gives  an  accurate  model  of  consolidation  of  saturated 
soils.  The  spring  represents  the  compressible  soil  skeleton,  the  fluid  represents 
the  pore  water  that  fills  the  soil  voids,  and  the  valve  represents  the  permeability 
of  the  soil.  It  is  easily  seen  that  the  rate  of  deformation  depends  upon  the 
soil's  permeability  (i.e.  how  much  the  'valve1  is  opened).  For  coarse-grained 
soils  (sands,  gravels)  the  permeability  is  so  high  that  the  deformation  response 
is  essentially  instantaneous.  Thus,  time-dependent  consolidation  effects  are 
generally  observed  only  in  fine-grained  soils  (silts,  clays).  The  entire  process 
can  be  made  precise  by  introducing  the  following  stress  concepts  (here,  in 
contrast  to  later  equations,  compression  is  taken  as  positive): 
t  :  the  total  stress  due  to  the  applied  load 
o':  that  portion  of  t  carried  by  the  soil  skeleton 
p  :  that  portion  of  t  carried  by  the  pore  water 

At  any  time  after  the  load  is  applied: 

t  =  o'  +  p  (1) 

Terzaghi  realized  that  the  deformation  of  the  soil  depended  directly  upon 
o'  and  not  t.  He  called  o'  the  effective  stress,  and  the  excess  pore  pressure 
p  the  neutral  stress  (since  it  does  not  directly  affect  the  soil's  deformation). 
This  concept  of  effective  stress  is  central  to  the  study  of  soil  mechanics,  whether 
consolidation  is  present  or  not. 

Terzaghi  developed  a  simple  model  [14]  for  consolidation  of  soil  layers, 
subject  to  the  following  assumptions,  see  Fig.  2. 

1.  The  consolidating  layer  is  horizontal,  of  infinite  extent  (laterally)  and 
of  constant  thickness  h. 
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2.  The  permeability  coefficient  (k)  and  volume  compressibility  (my)  are 
constant  throughout  space  and  time.  The  volume  compressibility  my 
represents  the  ratio  of  volumetric  strain  to  applied  effective  stress. 
Since  the  consolidating  layer  is  of  infinite  lateral  extent,  the  vertical 
strain  is  the  volumetric  strain. 

3.  The  pore  water  drains  only  in  the  vertical  (z)  direction. 

4.  The  time  rate  of  compression  depends  only  upon  low  soil  permeability: 
visco-elastic  properties  of  the  soil  skeleton  are  not  considered. 

5.  The  fluid  obeys  Darcy's  Law:  flow  is  proportional  to  the  gradient  of 
pore  water  pressure. 


v  =  -  r-  pM  (2] 

Tw 

6.  Strains  are  small  compared  to  unity. 

7.  The  applied  load  T  is  constant  for  all  time.  Thus,  since  T  is  given, 
if  the  pore  pressure  p  is  known,  so  is  the  effective  stress  o',  Eq.  (1). 

With  t  ese  assumptions,  Terzaghi  derived  a  differential  equation  for  this 
one  dimensional  case: 


p  =  vvp’ii  1=3  (31 

This  equation  is  identical  to  the  heat,  or  diffusion  equation,  and  it  can 
be  solved  using  separation  of  variables  for  various  boundary  conditions.  Fig.  3. 
Unfortunately,  many  of  the  assumptions  made  are  unrealistic  enough  to  warrant 
a  more  general  theory.  In  particular,  material  properties  are  not  constant,  and 
consolidating  layers  are  not  of  uniform  width  or  of  infinite  lateral  extent. 
However,  Terzaghi's  theory  has  been  widely  used  in  estimating  foundation  settle¬ 
ments  and  consolidation  rates. 
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hi  1941,  Biot  [4]  extended  Terzaghi's  consolidation  theory  to  the  general 
three-dimensional  case,  and  considered  loads  that  varied  with  time.  If  this 
extended  theory  is  restricted  to  fully  saturated  soils  where  the  applied  load  is 
constant  over  time,  the  generalized  model  is  governed  by  the  same  equation  as 
the  generalized  heat  (or  diffusion)  problem: 


P  = 


(sun  over  i) 


(4) 


As  in  the  one-dimensional  case,  this  equation  can  be  solved  using  separation 
of  variables  for  simple  geometries  and  boundary  conditions.  In  later  papers 
[5,15,16,17],  Biot  extended  his  consolidation  theory  to  include  effects  of 
anisotropy,  inhomogeneity,  and  more  general  boundary  conditions.  (Unfortunately, 
his  choice  of  notation  and  physical  constants  underwent  a  number  of  changes. 
For  simplicity  and  consistency  in  later  comparisons  with  finite  element  solutions, 
Biot's  results  are  paraphrased  in  the  following  development). 

The  resulting  system  of  differential  equations  for  elastic  consolidation  of 
an  anisotropic  soil  mass  can  be  summarized  as  follows: 

и.  :  soil  displacement  vector 

e.j  :  strain  tensor  (tension  positive) 

o'.j  :  effective  stress  tensor  (tension  positive) 

F-  :  body  force  vector  (per  unit  mass) 

к. .  :  permeability  coefficient  tensor 

^ijkl  5  elasticity  tensor  in  terms  of  effective  stress 
6..  :  Kronecker  delta 

p.  :  pore  water  pressure  (compression  positive) 

Vj  :  Darcy  velocity 


Strain-Displacement  Relations: 


eiJ  =  2<ui,j 


Equilibrium  Equations: 


(a. .  -  5. .  p) ,  +  p  F. 

i  j  *  j  J  s  1 


Effective  Stress-Strain  Relation: 


°ij  =  Ei  jki  ekil 


Darcy's  Law: 


r.  —  (k.  .  p, 

1  Yw  *1 


Equation  of  Continuity: 


e . .  «■  v . 

n  i,  j 


Biot  solved  these  equations  for  such  three  dimensional  cases  as  a 
rectangular  load  distribution,  Fig.  and  a  soil  with  an  impervious  top  layer. 
He  also  introduced  analytic  techniques  for  solution  of  a  variety  of  consolidation 
problems,  modelling  the  consolidating  layer  as  an  elastic,  semi-infinite  half  space. 
Although  the  mathematical  effort  required  >o  solve  these  equations  is  formidable, 
the  results  are  of  limited  pracu.  >•'  i  •  po: ; -u;  .  v.cee  assuming  infinite  depth 
for  a  soil  layer  can  lead  to  serious  overestimates  of  total  and  differential 


settlements.  The  utility  of  Biot's  ro, u 


.alytical  solutions,  but 


in  the  development  of  a  fairly  general  three  dimensional  theory  of  consolidation. 

3.2  Some  Finite  Element  Models  for  Consolidation 

One  of  the  main  reasons  for  the  widespread  use  of  finite  element  models 
in  mechanics  is  that  the  method  can  be  used  on  -oolern ■■  with  complex  geometries, 


variable  boundary  conditions,  and  non-hom ogeneous  material  properties.  Since 
these  difficulties  prevent  analytic  solution  of  Biot's  consolidation  equations  for 
many  useful  problems,  finite  element  procedures  have  become  an  Important  tool 
in  modelling  soil  consolidation  behavior. 

Sandhu  and  Wilson  [6]  published  one  of  the  first  finite  element  models 
for  soil  consolidation  in  1969.  This  model  combined  two  existing  finite  element 
models,  one  for  the  plane  strain  structural  problem,  and  the  second  for  the 
solution  of  the  two-dimensional  diffusion  equation,  see  Eq.  (9).  A  generalization 
of  Biot's  statement  of  Darcy's  Law  was  used,  which  included  the  effect  of  body 
forces  on  the  pore  water: 

p  :  water  density 

:  permeability  coefficient  tensor 

vi  *  kij(p»  j  *  vy s  0  (l0) 

Sandui  and  Wilson  used  a  variational  approach  for  the  derivation  of 
element  equations,  combining  finctionals  corresponding  to  the  two  coupled 
problems: 

0(,)  **ii  - 2p  *  **i . i  ♦*'  *vi  *"*j 

-  2ps  F,  *  Uj  ♦  g'  *  p,,  *  R,F,)dV  (11) 

t 

where  g'  =  1  and  a*b  =  J  a(t)  b(t-s)  ds 

o 

In  order  to  accomodate  traction  and  flow  boundary  conditions  (natural 
boundary  conditions  for  the  two  ooupled  problems),  the  following  terms  must  be 
added  to  the  functional  G(t)  to  obtain  the  desired  functional  F(t)  for  the  problem: 


Tj  :  components  of  prescribed  traction  vector 

Sj  :  portion  of  boundary  where  traction  is  prescribed 

Q  :  prescribed  normal  flux  on  boundary 

S2  •  portion  of  boundary  where  flux  is  prescribed 

n.  s  direction  cosines  of  outward  normal 

F(t)  =  G(t)  -  2  J  T.  *  u.  ds .  -  2  J  g'  *  Q  *  p  ds-  (12! 

s  1  1  1  s  * 

where  Tj  =  (a^.  -  6 ^  - p )n •  on  Q  •.  (n.  on  S2 

Sandhu  and  Wilson  evaluated  the  cortvokn.  integrals  using  a  simple  two 
point  forward  difference  formula  to  obtain  a  fully  explicit  time  marching  scheme. 
A  mixed  interpolation  model  was  used  on  the  displacement  and  pore  pressure 
unknowns:  displacements  were  interpolated  using  quadratic  shape  functions,  and 
linear  shape  functions  were  used  to  interpolate  the  pressure  unknowns.  Triangular 
elements  were  used  to  define  the  mesh:  there f  re  the  elements  incorporated 
a  six  node  linear  strain  triangle  for  the  structural  (displacement)  problem  and 
a  three  node  linear  triangle  for  the  fluid  proolern  (pore  pressure).  This  gives 
a  total  of  fifteen  degrees  of  freedom  for  each  element. 

Sandhu  and  Wilson  applied  the  finite  element  model  to  two  problems: 
Terzaghi's  one-dimensional  soil  column  and  Biot's  rectangular  strip  load  on  an 
elastic  half-space.  In  both  problems,  excellent  agreement  between  the  finite 
element  model  and  the  analytic  solutions  was  obtained,  see  Figs.  5  and  6.  Some 
discrepancy  can  be  seen  in  the  comparison  with  Bi-t's  c trip  load  solution,  but 
these  differences  can  be  attributed  to  tn«.  id,  t  tn.it  the  finite  element  mesh 
(like  the  soil  layers  being  mode lb  J)  is  of  limited  extent,  unlike  Biot's  elastic 


half  space. 
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A  later  work  by  Ghaboussi  and  Wilson  [S]  will  be  discussed  in  a  subsequent 
section. 

In  the  discussion  of  Biot's  generalization  of  Terzaghi's  consolidation  model, 
It  was  noted  that,  under  some  circumstances,  the  heat  or  diffusion  equation  can 
be  used  to  model  soil  consolidation.  Christian  and  Boehmer  [7]  developed  a 
displacement-based  finite  element  model  for  soil  consolidation,  and  compared 
results  for  the  finite  element  analyses  with  a  consolidation  analysis  based  on 
the  diffusion  equation.  Recall  that  the  diffusion  form  of  the  consolidation 


equation  (appropriate  when  the  applied  load  is  constant  over  time)  can  be  written: 


-v  F.jj 


The  volume  compressibility  m,  however,  is  not  a  material  constant.  It 
depends  on  the  type  of  analysis.  For  instance,  in  the  one-dimensional  case,  the 
volumetric  strain  is  simply  the  vertical  strain  e,  so  m  is  the  reciprocal  of 
the  constrained  elastic  modulus  (c  *  e  *  0).  However,  in  three  dimensions, 

*  J 

the  volumetric  strain  is  c  =  e  ♦  e  ♦  e_,  and  the  value  of  m  (and  hence  c  ) 

v  n  y  i  V  V 

should  be  modified  accordingly.  Christian  and  Boehmer  derived  correct  expres¬ 
sions  for  the  consolidation  coefficient  cy  and  total  stress  x  for  an  isotropic  soil 
mass  for  one,  two,  and  three-dimensional  cases  (see  Table  I).  With  these 


definitions  for  c  and  x,  Biot's  consolidation  theory  can  be  summarized  in  the 
general  equation: 

cv  p,jj  =  P  “  ^  (13) 

(Note  that  if  the  applied  load  is  constant  over  time,  x  =  0,  and  Eq.  (13)  becomes 
the  diffusion  equation). 

Christian  and  Boehmer  formulated  the  consolidation  problem  so  the  only 
unknowns  were  the  displacements  u.  Because  of  this  simplification,  the  near- 


Table  1 


Dimension 

(1) 


T 

(2) 


incompressibility  of  the  soil-water  system  causes  near-infinite  terms  to  occur 
in  the  equation  set  (see  Zienkewicz  [10]  for  an  explanation  of  this  phenomenon). 
In  order  to  remove  this  problem,  Christian  and  Boehmer  re-introduce  the 
continuity  condition  in  terms  of  pore  pressures  and  volumetric  strain  into  the 
equation  set. 

Christian  and  Boehmer  compared  their  finite  element  model  with  experi¬ 
mental  data  from  consolidation  tests  in  a  triaxial  load  apparatus.  Their  results 
agreed  well  with  experimental  and  analog  solutions.  One  surprising  type  of 
behavior  was  discovered  using  the  finite  element  model:  pore  pressures  locally 
increased  with  time  for  some  of  the  triaxial  tests.  Although  this  effect  was 
short-lived,  it  was  important  because  the  approximate  diffusion  equation  solution 
cannot  model  this  phenomenon.  Apparently  the  average  consolidation  can  increase 
(a  global  effect)  while  ome  parts  of  the  sample  experience  a  decrease  in 
effective  stress  (a  local  increase  in  pore  pressure). 

The  most  useful  results  of  this  study  by  Christian  and  Boehmer  were 
twofold:  first,  they  showed  that,  except  for  early  stages  of  consolidation, 

diffusion  solutions  can  be  utilized  when  the  problem  is  simple  enough  to  make 
such  a  solution  meaningful.  Second,  they  derived  correct  expressions  for  the 
appropriate  volume  compressibility  my,  depending  upon  the  geometry  of  the 
analysis.  These  values  of  my  insure  that,  if  a  diffusion  solution  is  used,  it  will 
be  one  that  is  appropriate  to  the  problem  being  analyzed. 

However,  the  finite  element  model  proposed  by  Christian  and  Boehmer, 
because  it  did  not  directly  model  the  pore  pressure,  required  some  manipulation 
in  order  to  handle  incompressibility  of  the  water-soil  system  and  certain  types 
of  flow  boundary  conditions. 

The  finite  element  model  proposed  in  1971  by  Yokoo  [9],  et  al.  is  practially 
identical  with  that  proposed  by  Sandhu  and  Wilson  [6]  in  1%9  (the  latter  paper 


was  discussed  earlier  in  this  review).  Yokoo's  work  was  completely  independent 
of  Sandhu  and  Wilson,  and  was  more  general  than  the  earlier  work.  In  addition, 
the  development  of  the  finite  element  equations  was  presented  in  a  more  lucid 
form  by  Yokoo,  who  carefully  considered  admissibility  and  boundary  conditions 
in  the  analysis.  Another  difference  between  the  twc  model'-  is  that  Yokoo 
evaluated  the  convolution  integrals  that  account  for  the  time  marching  by  using 
a  step-by-step  method  originally  derived  by  Zienkiewicz  and  Parekh  [19], 

Yokoo  considered  two  examples:  the  'classic'  one-dimensional  problem 
solvable  using  Terzaghi's  theory,  and  a  more  practical  axisymmetric  problem. 
In  the  one-dimensional  problem  Fig  7,  excellent  agreement  between  the  finite 
element  model  and  Terzaghi's  diffusion  solution  was  obtained. 

The  axisymmetric  problem  modelled  by  Yokoo  was  that  of  a  uniform  load 
on  a  circular  plate  loading  a  uniform  clay  layer.  The  clay  layer  exhibits  both 
structural  and  h  draulic  anisotropy.  In  addition,  the  applied  load  is  not  constant 
over  time.  This  problem  is  of  interest  because  it  is  a  useful  approximation  to 
a  common  foundation  problem  and  many  of  these  characteristics  (the  lot  al 
distribution  of  load,  anisotropy,  and  tune-dependent  loading)  violate  the  assump¬ 
tions  of  Terzaghi's  theory.  Yokoo's  paper  also  contains  excellent  graphical 
interpretations  of  the  evolution  of  pore  water  pressure,  and  of  the  displacement 
of  the  soil  layer. 

Perhaps  the  greatest  advantage  that  th*-  finite  element  method  has  over 
classical  (analytic)  solutions  for  -nsoh  i.ttioi  :s  that  irregular  geometry  causes 
no  serious  difficulties  in  a  finite  element  analysis.  This  is  particularly  important 
in  Geotechnical  Engineering,  because  soil  deposits  are  generally  irregular, 
nonhomogeneous  and  anisotropic. 

Desai  and  Saxena  [10]  ta«e  advan’  go  of  this  powerful  feature  of  the 
finite  element  method  and  model  cunsoiidatu  n  ol  a  layered  soil  deposit  underlying 


Figure  7.  One  Dimensional  Consolidation  From  Yokoo's  Model  [9] 


a  three  building  system.  The  actual  geometry  of  the  building  site  is  shown  in 
Fig.  8. 

Desai  and  Saxena's  analysis,  like  the  others  considered  earlier,  is  based 
on  Biot's  consolidation  theory.  Their  model  uses  a  mixed  interpolation  scheme, 
with  linear  interpolation  for  pore  pressure  and  quadratic  approximations  for 
displacement.  A  centered-difference  (Crank-Nicoison)  time  stepping  procedure 
is  used  to  numerically  integrate  the  convolution  integrals  that  account  for  the 
temporal  dependence  of  the  variational  terms.  A  useful  feature  of  their  model 
is  the  ability  to  vary  the  length  of  the  time  step:  since  consolidation  solutions 
exhibit  behavior  resembling  exponential  decay,  the  time  steps  can  be  lengthened 
as  the  solution  asymptotically  approaches  the  steady-state  solution.  (The  parti¬ 
cular  problem  is  modelled  over  a  period  of  almost  twelve  years,  with  time  steps 
ranging  from  one  to  forty  days.) 

Because  soil  deposits  are  heterogeneous,  the  determination  in  geotechnical 
engineering  of  material  parameters  often  becomes  a  statistical  exercize.  For 
this  reason,  Desai  and  Saxena  analyzed  several  models  of  the  three-building 
foundation  problem,  each  with  slightly  different  material  properties  and/or  degree 
of  anisotropy.  The  results  they  obtain  are  very  interesting.  Fig.  9. 

Five  of  the  six  finite  element  analyses  can  be  se-m  to  closely  approximate 
the  actual  measured  consolidation  settlement  (shown  for  building  2).  The  one 
analysis  that  gives  a  poor  result  has  an  unlikely  type  of  anisotropy,  where  the 
vertical  permeability  exceeds  the  horizontal  permeability  by  a  factor  of  one 
hundred  (t^e  horizontal  permeability  is  generally  larger).  Desai  and  Saxena  also 
calculate  an  estimate  for  the  probable  ultimate  settlement  of  this  building,  using 
Terzaghi's  one-dimensional  consolidation  theory.  Although  this  theory  clearly 
does  not  apply  here,  it  is  extensively  used  in  similar  cases  to  obtain  settlement 
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Figure  8. 1  Geometry  of  Building  Site  Modelled  Q(Q 


estimates:  these  estimates,  in  general,  are  very  conservative.  However,  in  this 
case,  actual  settlements  are  250%  of  those  predicted  by  Terzaghi's  theory! 

Desai  and  Saxena's  work  is  characterized  by  a  practical  emphasis,  both 
in  terms  of  the  useful  problems  solved,  and  the  rules  of  thumb  proposed  for 
future  analyses.  With  this  work,  the  advantages  of  a  finite  element  analysis 
for  consolidation  are  more  fully  realized. 

Three  of  the  assumptions  inherent  in  Terzaghi's  original  consolidation 
theory  are  sufficiently  unrealistic  to  warrant  development  of  a  better  model: 

A.  The  soil  skeleton  is  a  linear  elastic  material  (i.e.  the  volume  compressi¬ 
bility  my  is  constant) 

B.  Vertical  Flow  (one-dimensional  behavior) 

C.  Applied  (total)  stresses  are  constant  over  time. 

The  various  finite  element  models  discussed,  based  on  Biot's  consolidation 
theory,  have  removed  the  restrictions  due  to  the  second  and  third  assumptions. 
However,  none  of  them  have  attempted  to  model  the  soil  skeleton  as  an  inelastic 
material.  Thus  the  area  of  consolidation-related  research  of  greatest  current 
interest  involves  the  inclusion  of  inelastic  soil  effects.  In  addition  to  the  work 
reported  in  the  next  section,  the  reader  is  referred  to  [11,  12,  and  20]. 

SUMMARY  OF  GOVERNING  EQUATIONS  USED  IN  THIS  WORK 

For  the  sake  of  completeness  many  of  the  concepts  and  equations  given 
in  the  previous  section  are  repeated  here.  With  only  minor  changes  in  notation, 
sign  convention  and  theory  the  following  is  taken  from  reference  [8] ;  the 
significant  difference  is  the  use  of  bounding  surface  plasticity  theory  to  model 
soil  behavior. 

Throughout  the  following  sections  the  usual  convention  is  used  that  free 


indices  can  vary  over  their  ranges  and  repeated  indices  must  he  summed  over 
their  ranges;  commas  denote  differentiation. 


4.1  Concept  of  Stress 


The  pore  water  pressure,  total  (phenomological)  stress,  and  effective  stress 
are  denoted  respectively  as  h,  t and  o’...  Pore  water  pressure  is  taken  to 
be  positive  in  compression  while  the  mechanics  sign  convention  of  tensile  normal 
stresses  being  positive  is  used  for  tjj  and  aVj.  For  the  purposes  of  the  theoretical 
development  h  is  taken  to  be  the  total  pore  water  pressure  (i.e.,  including  the 
hydrostatic  pressure).  Later  in  the  discussion  of  the  finite  element  programs, 
means  for  treating  it  either  as  total  or  "excess'1  pressure  are  discussed,  the 
effective  stress  o'^  is  that  portion  of  the  total  stress  carried  by  the  soil  skeleton. 

The  relation  between  these  three  stress  quantities  is  [8] : 


(14) 


Not  all  authors  include  the  factor  a  in  the  above  equation.  It  appears  that  its 
presence  permits  consideration  of  the  fact  that  the  average  stress  contribution, 
over  a  unit  cell  of  soil,  due  to  the  pore  pore  water  pressure  may  be  less  than 
the  actual  pressure.  For  example,  consider  the  idealized  case  illustrated  in  Fig. 
10.  However,  for  actual  cohesive  soils  there  is  very  little  stress  transfer  through 
particle  contact,  thus  in  practise  o  should  be  very  nearly  unity.  In  fact  most 
researchers  [6,9,10]  do  not  include  the  quantity  and  Ghaboussi,  et  al.  [8]  set 
it  equal  to  unity  for  all  examples  considered.  Finally  if  a  is  included  in  this 
expression  and  one  wishes,  for  computational  purposes,  to  express  the  governing 
equations  in  variational  form  it  must  also  be  included  in  the  conservation  of 
mass  equation.  Its  physical  significance  in  this  second  equation  is  not  clear  and 
its  inclusion  would  appear  to  be  arbitrary.  Thus  equation  (14)  is  rewritten  with 


4.2  Conservation  of  Mass 


Denote  the  total  volume  of  water  that  has  flowed  out  of  a  unit  volume 

of  soil  by  V  ,  the  strain  of  the  soil  mass  by  e..  (tensile  strain  positive)  and  the 
w  IJ 

total  displacement  of  the  soil  by  u^  For  small  strains 


ij 


UJ.I» 


(16) 


It  is  assumed  that  the  soil  particles  and  the  pore  water  experience  an  elastic 
decrease  in  volume  due  to  an  increase  in  pore  water  pressure;  denoting  the 
corresponding  bulk  modulus  by  T ,  the  resulting  rate  of  volume  change  ish/T. 
(Alternatively  the  parameter  T  can  be  thought  of  as  a  "penalty  number"  used 
to  approximately  imp  se  an  incompressibility  constraint  on  the  water  and  soil 
particles.)  This  expression  embeds  the  assumption  that  the  mean  stress  component 
of  the  effective  stress  produces  no  significant  volume  change  of  the  soil  particles. 

For  a  completely  saturated  system  the  rate  of  volume  change  of  the  soil 
£..)  most  be  balanced  by  the  water  flowing  out  and  the  rate  of  volume  change 
of  the  water  and  soil  particles,  i.e. 


(17) 


In  [8]  the  factor  a  of  equation  (14)  multiplies  the  e~  term  in  eq.  (17).  The 
presence  of  a  in  eq.  (17)  is  necessary  if  it  appears  in  eq.  (14)  and  it  is  desired 
to  represent  the  problem  by  a  variational  statement;  however,  its  physical 
meaning  in  eq.  (17)  is  unclear. 


4.3  Water  Flow 


An  average  displacement  w.  of  the  fluid  relative  to  the  soil  is  defined 
such  that  WjnjdA  is  the  total  amount  of  fluid  that  crosses  the  area  dA  (outwa<  d 
normal  n.).  Thus  the  average  (Darcy)  velocity  oi  the  fluid  is 


(18) 


Of  course,  the  actual  velocity  u.  the  pores  of  t:  >  soil  is  much  higher.  Denote 
the  effective  permeability  tensor  of  tor  so:(,  1  d  density  and  component 

of  gravity  in  the  i  direction  as  kV,  p^  o:.  .hy.  The  term  effective 

permeability,  as  used  here,  is  equal  to  the  rat-  ;  permeability  coefficient 

commonly  used  in  civil  engineering  literal  i"-  ied  by  the  unit  weight  of 

water,  which  in  turn  is  equal  to  the  geo  i>  ;  r •- meability  used  in  physics 
divided  by  the  viscosity  of  the  water.  In  terms  of  this  notation,  Darcy's  law- 
governing  the  water  flcxv  is  expressed  as 


w.  =  v.  = 

i  i 


k*. . (h, . 
‘1  ) 


(19) 


For  this  application  it  is  assumed  that  toe  eii -  o  cx-r  meability  tensor  is  a 
constant  (i.e.,  its  components  are  not  strain  depen  a--- it  >  and  symmetric. 


4.4  Equilibrium 

Denoting  the  total  mate,  la:  density  ;i,  soil)  by  p,  the  linear 

equilibrium  equations  take  on  then  us...:  • 

T.  .  .  +  pg.  =  0  (20! 

4.5  Strain  -  Effective  Stress  Relations 

Soils  are  in  general  nonlinear,  meias;.  once,  the  stress  strain 

law  will  involve  some  type  of  hereditary  The  implementation  of 


such  a  relationship  usually  requires  a  step-by-step  (incremental)  solution.  For 
a  given  time  step  N,  this  relationship  can  be  expressed  in  the  form  (for 
convenience  the  increment  number  N  will  not  be  displayed  but  merely  implied). 


Ac' . .  =  D- 
i)  i 


ikt  AEkt  ♦  4a'i), 


(21) 


For  actual  use  in  a  finite  element  analysis  this  equation  would  be  written  in 

matrix  form  [2,21]  however,  for  the  purposes  of  equation  development  the  tensor 

notation  is  preferable.  In  general  the  tensor  of  incremental  properties  D..^  is 

a  function  of  the  solution  (i.e.,  Aa'~  and  Ae^),  and  thus  some  nonlinear  solution 

scheme  employing  iteration  is  usually  necessary.  In  many  cases  the  term  Ao'.. 

% 

is  dependent  only  upon  the  past  history  (0  -  t^j)  of  the  solution,  however,  in 
general  it  may  also  depend  upon  Ao'-j  and  Ae^.  For  the  current  bounding 
surface  application  it  is  zero.  Equation  (20)  is  sufficiently  general  to  accomodate 
the  model  of  interest  in  this  study,  the  bounding  surface  model  for  cohesive 
soil,  as  well  as  linear  elasticity  and  most  other  standard  and  advanced  models. 
For  the  actual  functional  form  of  for  bounding  surface  theory  the  reader 

is  referred  to  ref.  [2,22]. 


4.6  Boundary  Conditions 

For  the  sake  of  brevity  only  simple  boundary  condition  are  considered  at 
this  time  (i.e.  spring  and  convection  type  conditions  are  excluded).  At  every 
point  on  the  boundary  of  the  soil  mass  the  rates  of  either  the  traction  or 
displacement  components  will  be  given,  i.e. 


tj  =  Tj.  n.  or  u.  given 


(22) 


28 


The  components  of  the  unit  normal  to  the  surface  are  n..  In  addition  the  rate 
of  the  pore  water  pressure  "h"  or  the  total  flow  of  water  (Q  =  wm.)  will  be 
known,  he. 

h  or  Q  g  ven  (23) 


4.7  Incremental  Equations 

Because  of  the  time  dependent  nature  of  the  problem  a  step-by-step 
solution  scheme  is  required.  Thus  the  variables  u  ->xpressed  in  an  incremental 
form,  i.e.  for  time  N  u.^  =  ,  ♦  An.,  :  v.ously  noted  the  subscripts 

N  on  the  incremental  values  not  <  .  .played.  It  is  now  necessary 

to  express  equations  (15)  -  (21)  in  inc.i e-aeota! 

Consider  first  eq.  (17).  Recall  uk.,  ,  o  it;  he  total  volume  of  water 
that  has  flowed  out  of  a  unit  volume  of  soil;  to.*  ... *c  of  this  quantity  is  simply 
related  to  the  average  fluid  velocity,  i.e. 


w. 

1 , 1 


Substituting  this  expression  into  eq.  (17) 

•  •  h 

Wi,i  =  "  eit  ‘  r 

Integrating  the  above  equat  -  .'m- 

constant  then  some  form  of  .i,  rv<.  ’ 
last  term): 


(24) 


(25) 

■:  yields  (if  T  is  not  a 
■Id  be  required  for  the 


Aw. 


i .  i 


(26) 

iting  a  weighted  residual 


(If  one  prefers  the  above  step  can  be  ;-e 


Integrating  eq.  (19)  over  the  time  interval  gives: 


A*i  =  -  /  k*^  (h,.  -  pfg.)  dt  (27) 

Vi  3  3 

It  is  assumed  that  k*„  is  a  constant  (permitting  it  to  be  state  dependent,  hence 
implicitly  time  dependent,  would  only  slightly  complicate  the  analysis),  i.e. 


In  order  to  accomodate  possible  centrifuge  applications  the  gravity  term  is 
considered  to  be  time  dependent.  The  two  integrals  on  the  right  are  now 
approximated  using  numerical  integration.  Trapezoidal  integration  is  used  on 
the  second  term,  while  the  more  general  rule 


/  F(t)dt  *  [(1-0)  Fn_1  +  SFjjAt  =  [Fn1  +  QAFlAt  (29) 
fN-i 

is  used  for  the  first  term.  Thus,  eq.  (28)  yields  (p{g.  =  ^  Kpfgj)N_!  +  (pfgj)N)): 

Awj  =  -  k*j  j  ^  +  9Ah,.  -  p{gj]At  (30) 

Values  of  6  of  0,  1  and  1/2  give  forward  integration,  backwards  integration  and 
the  trapezoidal  rule  respectively;  alternatively  if  one  prefers  to  discretize  time 
by  approximating  the  time  derivative  in  eq.  (19)  using  a  finite  difference  operator 
the  cited  values  of  9  correspond  to  using  a  forward  difference  (Euler's  method), 
a  backward  difference  and  a  central  difference  (Crank-Nicolson  or  mid-point 
method);  finally  if  one  prefers  the  weighted  residual  interpretation  [18]  these 
values  of  0  correspond  to  a  delta  function  weight  at  the  backward  point,  a  delta 
function  weight  at  the  forward  point  and  a  uniform  weight  respectively.  Values 


30 


of  0  <_  1/2  give  schemes  which  are  only  conc'itu  nally  stable  and  should  be 
avoided.  Theoretically  0-1/2  should  lead  rn  the  greatest  accuracy,  but 
practically  may  lead  to  oscillation  problems.  Ti  t  retically,  a  value  of  0  =  1 
eliminates  all  oscillate  bur  1  ,.o  cn  v  . ,  i  s-.ar?  ct  eristics.  HjenKiewicz 

[18]  suggests  a  compm  use  •>  -  cf  •  *1  alls  "Galerkin’s"  method  as 

it  bears  a  resemblance  to  Gaierkm's  ■  ghted  rev  j,  method  for  boundary  value 
problems.  The  unlimited  possioi  uies  offered  hy  r,,  ccr -order  integration  formulas 
(e.g.,  improved  trapezoidal  etc  )  are  ->ot  expt  >re  :  here. 

The  incremental  forms  of  eqns.  ii  >).  mu  v20)  are  found  imply  by 

writing  the  equations  at  t^  ^  and  t„,  im  sub  ,  .g.  Equations  (22)  and  (23) 
are  converted  to  incremental  fore  t. .  r  the  interval.  The  results 


along  with  eqs.  (21)  (26)  and  (30)  are  sum.i 

combined): 

Field  equations: 

low  (eqs.  (15)  and  (21)  are 

Axi)  =  DI  jkl  Atk?.  •‘  !i  \  ; 

-  i 

(31) 

Av .  =  -  k  v  .  [h. ,  .  *  \i  , 

i  i j  N- 1 , 

'  "  i "  i  J 

(32) 

Ae . .  =  i( Au .  .  +  Au  ) 

» J  2  >,)  J  > 1 

(33) 

At..  .  +  A(pg.)  0 

1  )  f  1  ) 

(34) 

-  “■»  *r 

boundary  conditions: 

(35) 

M 

o 

-i 

r> 

o 

— i 

(36) 

and 


IQ  or  Ah  given 


(37) 


This  set  of  equations  con  be  reduced  in  number  by  expressing  them  in 
terms  of  primary  dependent  variables  Au^  and  Ah.  Equation  (33)  is  substituted 
into  eq.  (31)  and  the  results  into  eq.  (34),  and  eqs.  (32)  ^nd  (33)  are  substituted 
into  eq.  (35)  to  yield: 


<Dijki  tuk,i),i  -  sij  *  to'ij0j.  *  - 0 


(kVhN-i,j  ♦  9ah-j  -  OfSjU.i At  -  4ui,i  -  r  = 0  <»> 

The  corresponding  boundary  conditions  are  found  by  substituting  eqs.  (31),  (32) 
and  (33)  into  eqs.  (36)  and  (37): 

ATj  *  <Dijk!  AuM  ‘  ^  {ij  *  to’iJ  )ni  °r  Auj  «iven  (,,0> 


+  9Ah,.  -  Pfgjln^  At  or  Ah  given  (41) 

Thus  eqs.  (38)  and  (39)  are  the  final  form  of  the  governing  incremental 
equations  (for  time  step  N),  subject  to  the  boundary  conditions  of  eqs.  (40)  and 


4.8  Variational  Statement  of  the  Problem 

The  solution  of  the  boundary  value  problem  given  by  eqs.  (38)  -  (41)  will, 
except  in  the  simplest  cases,  require  numerical  analysis.  A  finite  element 
solution  of  these  equations  can  be  formulated  directly  by  applying  Galerkin's 


c-.r 


weighted  residual  met  nod  or  .iirernatively  by  seeding  a  stationary  point  for  ar. 
appropriate  variational  staterru  of  the  problem;  me  latter  method  is  used. 

Using  variational  calculus  concepts  it  is  a  s.  i.pie  matter  to  construe;  a 
variational  statement  that  is  equivalent  to  these  equations;  it  has  the  following 
form: 

6F  =  0  (42) 


where 


I  (Ah)2 


-  9  --  k*.  Ah, .  tU  r\' 

^  i)  )  .  i i ; 

-  4t  -  h'8.  4ui|dv 

j  AT-  Au-  dS  f  As<  Ah 

c  l  i  , 


(43) 


Tite  volume  of  me  soil  mass  is  d«-n. iy;  \  md  the  surface  areas  over 
which  AT-  and  AQ  are  specified,  b>  Sj  and  S,.  *  , a  alternative  course  of  action 
is  to  bypass  the  incremental  differential  eq..  ,*.,n  .a.  1  to  obtain  eq.  (43)  directly 
from  a  variational  statement  of  the  type  ;-  •>  r  e  q.  (12)  and  a  numerical 
approximation  to  the  convolution  integrals;  tin  ,q»  ■  .  would,  however,  require 
a  somewhat  different  numeric,  tte  ’  '-ei..  •  tding  surface  plasticity 

model.) 

5.  FINITE  ELEMENT  ANAL  A  ' 

For  the  two-  and  three -dm  >  •!  analyses  developed  as 

part  of  this  project,  standard  i-  eametric  elements  are 


used.  This  selection  was  based  on  a  consideration  of  the  ease  of  data  preparation 
and  to  a  lesser  extent  on  the  results  of  an  informal  report  by  Professor  Segerlind 
of  Michigan  State  which  indicates  that  for  time  dependent  problems  the  low 
order  elements  experience  fewer  oscillation  problems  them  the  higher-order  ones. 
Because  the  steps  required  to  proceed  from  eq.  (42)  to  the  finite  element 
equations  are  well  documented  (e.g.,  see  [18])  only  a  few  special  considerations 
are  discussed  here. 

For  a  linear  elastic  material  is  constant  and  the  analysis  proceeds 

simply  in  a  step-by-step  fashion.  However,  for  a  cohesive  soil  characterized  by 
the  bounding  surface  plasticity  model,  this  tensor  is  highly  dependent  upon  the 
solution  and  hence  each  incremental  analysis  is  decidedly  nonlinear.  The 
approximate  Newton-Raphson  method  used  to  solve  the  nonlinear  problem  is 
throughly  discussed  in  [23].  The  two  programs  are  written  so  that  a  user 
specified  value  of  B  between  0  and  .5  is  used  to  select  an  approximate 
Newton-Raphson  method  on  the  spectrum  from  "tangent  stiffness"  to  "successive 
approximations"  [23,24,25]. 

The  one  characteristic  of  the  problem  that  requires  some  care  is  the 
handling  of  the  near-incompressibility  of  the  soil  when  it  is  in  a  saturated 
condition.  The  most  general  procedure  for  avoiding  the  accuracy  and  round-off 
error  problems  associated  with  the  finite  element  analysis  of  nearly  incompressible 
materials  is  the  use  of  a  "mixed  formulation"  analysis  [19,26].  Its  use  is  natural 
for  soil  consolidation  problems  because  the  additional  mean  pressure  variable, 
needed  in  the  mixed  formulation,  is  already  included  in  order  to  describe  the 
flow  problem,  i.e.,  eq.  (42)  is  a  natural  "mixed"  statement  of  the  problem. 

The  use  of  the  mixed  formulation  for  isoparametric  elements  must, 
however,  be  done  with  considerable  care.  The  problem  is,  if  the  near¬ 
incompressibility  condition  is  applied  point-wise,  the  elements  "lock-up",  i.e., 


become  rigid  and  will  not  deform  [181.  To  avoid  this  problem,  nmtr- 
incompressibility  condition  must  not  be  satisfied  point-wise  but  onlv  m  souse 
average  sense.  For  a  first  order  element  only  the  average  volume  change  lor 
tlie  element  can  be  made  zero.  This  is  accomplished  if  the  term  in  eq.  ('ill 
which  measures  the  volume  change,  measures  only  average  volume  change  (over 
the  element)  and  not  point-wise  change;  the  term  in  question  is  Ah  Am..  In 
order  to  measure  only  average  volume  change,  the  element  approximation  for 
Ah  in  this  term  must  be  a  constant.  However,  the  admissibility  condition  lor 
Ah  which  arises  out  of  the  presence  of  the  term  Ah,-  Ah,j  requires  that  its 
approximation,  for  this  second  term,  be  continuous  between  elements.  This 
incongruency  can  be  easily  dealt  with  by  using  a  two  field  approximation  for 
Ah.  The  first  (used  in  the  term  Ah  Ac--)  is  constant  for  each  element  and 
thus  not  continuous  across  element  boundaries,  the  second  is  continuous  and  is 
used  for  all  other  terms  in  eq.  (43).  The  two  fields  are  related  to  a  common 
set  of  nodal  unknowns.  The  continuous  field  is  defined  by  the  node  point  values 
and  first  order,  isoparametric  shape  functions.  The  constant  element  value  for 
the  discontinuous  field  is  defined  to  be  the  average  of  the  values  lor  the  nodes 
describing  the  element,  thus 


fe^i  *i 


Ah(2)  =  AH-  N. 


(45) 


Where  ND  is  the  number  of  nodes  defining  the  element  (4  or  Is),  N.  are  the 
first  order  isoparametric  shape  functions,  AH.  are  the  node  point  values  of  An 
for  the  nodes  defining  the  element  and  I.  =  1.  When  me  prefers  to  us.*  'h<* 
continuous  approximation  for  Ah  in  ail  terms  then  I.  is  replaced  by  i  ei 

the  programs  written  to  evaluate  the  analysis,  the  user  can  choose  between 
these  two  alternatives  by  means  of  a  simple  input  code. 
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The  isoparametric-Laplacian  grid  generation  scheme  given  in  [27]  has 
been  modified  in  order  to  replace  the  iterative  solution  by  a  direct  solution; 
the  reduction  in  computational  cost  for  this  step  in  the  analysis  is  dramatic. 
For  the  three-dimensional  code,  this  grid  generation  scheme  has  been  generalized 
to  produce  meshes  consisting  of  8-node  brick  elements. 

For  unsaturated  or  "ideal  drained"  conditions  the  h  variable  is  dropped  to 
reduce  the  number  of  unknown's  per  node  by  one.  However,  for  "ideal  undrained" 
conditions  and  a  saturated  soil,  the  h  variable  is  retained  to  facilitate  modeling 
the  resulting  near-incompressibility  (by  means  of  the  mixed  formulation  as 
explained  above). 

The  implementation  of  the  bounding  surface  model  followed  directly  the 
instructions  given  in  [22].  Reference  [22]  is  a  revision  of  a  portion  of  reference 
[2]  to  reflect  a  complete  recoding  of  subroutine  CLAY  and  its  affiliated 
subroutines.  CLAY  was  recoded  for  the  sake  of  clarity,  to  incorporate  some 
recent  minor  changes  in  the  model,  to  improve  numerical  efficiency  and  to  take 
advantage  of  the  structured  programming  concepts  of  "Fortran  77."  For  the 
sake  of  illustrating  the  implementation  procedure  (according  to  the  well- 
documented  instructions  in  [22])  of  the  bounding  surface  model  in  finite  element 
codes,  no  changes  whatsoever  were  made  in  CLAY  and  its  subroutines  for  this 
application.  As  a  result  one  unused  subroutine  is  retained  and  there  is  duplicate 
input  for  the  quantity  T.  In  order  to  simulate  the  modification  of  existing 
programs,  the  two  finite  element  programs  developed  for  this  project  were  first 
written  as  relatively  general  incremental-iterative  nonlinear  programs  and  then 
the  CLAY  subroutine  was  included  (by  means  of  two  simple  call  statements  per 
program)  as  a  modular  unit. 

The  notation  used  in  the  following  discussion  of  this  implementation  is 
the  same  as  used  in  [22].  Subroutine  RPROP  is  called  from  subroutine  PROPTY 


and  reads  the  parameters  which  describe  the  bounding  surface  model,  .'or 
convenience  the  combined  bulk  modulus  of  the  soil  particles  and  pore  water  r 
(not  really  a  parameter  of  the  plasticity  model)  is  read  and  stored  separately 
by  PROPTY;  the  read  for  T  in  RPROP  is  a  duplication  and  is  not  necessary. 
The  values  of  void  ratio  eQ  and  preconsolidation  pressure  pQ  in  array  STOR  are 
initialized  for  each  element  in  Subroutine  GEOM.  The  STOR  array  for  each 
element  is  included  in  the  "BLK6"  records  stored  on  unit  2. 

The  analysis  is  a  mixed  formulation  (see  Section  5)  and  thus  in  the  CALL 
to  CLAY  (from  PROPTY)  KIND  =  0.  The  finite  element  program  supplies  h  , 
and  Ah^  to  CLAY  through  the  "CALL".  Because  in  this  application  it  was 
found  convient  to  store  T  in  the  main  program,  the  quantity  GAM  in  the  CALL 
is  not  used.  The  combining  of  the  arrays  [Dl^j  K  j  and  [D]  N  K  ;  (according 
to  eqs.  (17)  and  (18)  of  [22])  is  done  in  PROPTY  immediately  after  the  CALL 
to  CLAY.  The  reversing  of  the  sign  convention  for  the  normal  stresses  and 
strains  (and  for  the  2-D  program,  the  expanding  of  the  two-dimensional  stress 
and  strain  vectors  to  three-dimensional  form)  is  done  just  prior  to  the  CALL 
to  CLAY.  For  this  small  deformation  analysis  LARGE  =  0.  LOCIT  is  set  equal 
to  ITMAX  used  for  the  global  iteration  and  ERMAX  is  set  equal  to  10  times 
the  value  used  in  the  global  iteration.  TH1  has  the  value  of  .*>  as  used  in  the 
main  program.  In  all  cases  IDIM  =  3  (this  is  true  for  plane  strain  conditions 
as  the  finite  element  analysis  calculates  oz).  Information  concerning  use  of  the 
programs  is  to  be  found  in  [28]  and  [29]. 

Tw  o-dimcreionaJ  element  matrices: 

For  plane  strain  conditions  to  exist,  the  only  off  diagonal  term  in  the 

k*..  tensor  that  can  be  non-zero  is  k*  (k*  ); 
ij  xy  ** 


this  coefficient  is  denoted  as 


Expressing  the  displacement  approximations  in  terms  of  their  node  point 
values  and  the  shape  functions,  using  eqs.  (44)  and  (45)  in  the  appropriate  terms 
and  differentiating  with  respect  to  the  node  point  unknowns  yields  the  element 
matrices.  The  terms  arising  because  of  the  presence  of  the  Ah  variable  are 
explicitly  given  below: 

52Q7  =  //{(  JAUj  ♦  I  ]AVj  +  [-  ^  (F-  +  pNi )] AHj  ♦  [  ]}RdA 


32r=JJ{t  1AU-  +  [  ]AV.  +  [-  /GjJAH.  +  [  ]}RdA 

i  A  ‘  1  1 

aic’tfh-r  (f,  +  pNjijAUj  ♦  i-^Cjiw, 

i  A 

+  {-  T  Ni  Nj  -  0  At[k*ll  FiFj  +  k*12(FiGj  +  FjGi)  +  k*22  GjCj] }AHj 

-  Attfk#ll(^-lk  Fk  -  Pf«l>  *  k#12(^-lk  S  -  Pf82>lFi 

+  [k*12(HS[_ik  Fk  -  pfgj)  +  k*22%-Ik  Gk  "  pfS2)lc,il  RdA 

The  terms  not  shown  are  identical  to  those  for  conventional  stress  analysis.  For 
plane  strain  conditions  R  =  1,  p  =  0,  while  for  axisymmetry  R  =  r  and  p  =  -. 
The  x(r)  and  y(z)  derivatives  of  the  shape  functions  (N^)  are  denoted  respectively 
by  F.  and  G.. 

i  i 

Three-dimensional  element  matrices: 

The  terms  in  the  element  matrices  arising  from  the  presence  of  the  Ah 
variable  are  given  below:  (Note  that  AW^  is  the  change  in  displacement  in  the 
z-direction,  not  the  average  fluid  displacement  W-  Also,  the  derivatives  of  the 
shape  functions  with  respect  to  the  global  coordinate  directions  x,  y  and  z  arc 


The  empty  brackets  indicate  terms  which  are  identical  to  those  found  in  a 
conventional  stress  analysis. 

6.  EXAMPLES 

During  the  check-out  phase  of  the  code  development  numerous  example 
problems  were  analyzed.  Results  from  two  of  these  analyses  are  given  in 
Figures  11  and  12. 

The  first  example  is  a  generalization,  to  include  the  compressibility  of 
the  soil  particles  and  pore  water  (T<®  in  eq.  (17)),  of  the  one-dimensional 
Terzaghi  problem.  In  Figure  11  the  finite  element  predictions  for  the  deflection 
at  the  surface  are  compared  to  the  exact  results  taken  from  [8] .  It  should  be 
noted  that  the  interpretation  of  the  results  given  in  [8]  must  be  done  with 
care.  The  contents  of  Figure  1  give  the  impression  that  the  solution  only 
depends  on  the  parameter  M/Cy,  whereas  it  is  easy  to  show  that  it  depends 
instead  on  the  quantity  (the  terms  are  defined  in  [8] ).  The  results  given 
in  Figure  1  of  [8]  appear  to  have  been  run  for  the  case  of  k/n  =  1.0  and  thus 
the  ambiguity  caused  no  problem. 

The  second  example  considered  the  uniform  loading  of  a  soil  layer  that 
is  free  to  drain  both  at  the  surface  and  into  a  central  sand  drain.  The  finite 
element  mesh  used  in  the  analysis  is  illustrated  in  [28].  Figure  12  compares 
the  predictions  for  the  surface  displacement,  at  a  radius  20  times  that  of  the 
drain,  to  the  results  given  in  [8], 
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